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Abstract 

The recent liberalization of the electricity and gas markets has resulted in the growth of 
energy exchanges and modelling problems. In this paper, we modelize jointly gas and elec- 
tricity spot prices using a mean-reverting model which fits the correlations structures for the 
two commodities. The dynamics are based on Ornstein processes with parameterized diffusion 
coefficients. Moreover, using the empirical distributions of the spot prices, we derive a class of 
such parameterized diffusions which captures the most salient statistical properties: stationarity, 
spikes and heavy-tailed distributions. 

The associated calibration procedure is based on standard and efficient statistical tools. We 
calibrate the model on French market for electricity and on UK market for gas, and then we 
simulate some trajectories which reproduce well the observed prices behavior. Finally, we illus- 
trate the importance of the correlation structure and of the presence of spikes by measuring the 
risk on a power plant portfolio. 
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, 1 Introduction 

■ The recent deregulation of energy markets has led to the development in several countries of market 
places for energy exchanges. Consequently, understanding and modelling the behavior of energy 

£f~^ | market is necessary for developing a risk management framework as well as pricing of options. 
CN ■ Many derivatives on both electricity and gas spot (and futures) prices are traded. Understanding 
the correlation structure between both energies is a significant challenge. For instance, spark spread 



q 

options are commonly traded in energy markets as a way to hedge price differences between elec- 
tricity and gas prices or are used in order to price projects in energy (see [12] for an introduction). 
Thus, modelling jointly the evolution of gas and electricity prices is a relevant issue. 

Numerous diffusion-type and econometric models have been proposed for electricity and gas 
spot prices. In energy markets, spot price dynamics are commonly based on Ornstein processes, 
which are the classical way to model mean-reversion. Geometric models represent the logarithmic 
prices by a sum of Ornstein processes with different speeds of mean reversion whereas arithmetic 
models represent the price itself (see for instance [21] for a geometric model). Also, equilibrium 
models ([2] and [T5]) have been investigated in order to reproduce price formation as a balance 
between supply and demand. The main drawback of such model is that they do not reproduce the 
autocorrelation structure of a commodity and the cross-correlation structure between commodities. 
In [13], a markov jump diffusion is investigated for electricity spot prices. Though, it properly 
represents the spiky behaviour of spot electricity prices, the process reverts to a deterministic mean 
level whereas it usually reverts to the pre-spike value on data. Moreover applied to electricity and 
gas spot prices, it does not capture the autocorrelation and cross-correlation struture observed on 
data. 
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Another class of spot price dynamics is represented by multifactor models. Several authors (see 
|14j . [1], [9], [20], [22] among others) have investigated this kind of diffusion. The logarithmic prices 
or the price itself is represented by a sum of Ornstein processes in order to incorporate a mixture 
of jump variations and "normal" variations. For instance, in [20] the deseasonalized spot price or 
log-spot price X(t) is given by: 

X(t) = Y 1 (t) + Y 2 (t) 

where 

dY i (t) = -\ i Y i (t)dt + dL i (t), i = 1,2. 

The Ornstein Uhlenbeck (OU) component Y\ is responsible for the normal variation and is assumed 
to be Gaussian, i.e. L\(t) is a Brownian motion, whereas Y 2 is the Levy driven OU component 
responsible for spikes, i.e. L 2 {t) * s a jump Levy process. In this kind of framework, the difficulty 
is to detect and filter the spikes in order to estimate the jump part. Several methods have been 
proposed to circumvent this problem (see e.g. |20j and [3]). For instance, |14] presents a similar 
model for the spot price process that is the exponential of the sum of an Ornstein-Uhlenbeck and 
an independent mean reverting pure jump process, derives its associated forward curves and finally 
proposes to calibrate it to the observed forward curve at time t = 0. In order to calculate premia 
of call and put options as well as path-dependent options several approximations to the probability 
density function of the logarithm of the spot price process at maturity are done. In [BJ, the following 
spot price dynamics for two energies A and B are proposed 

m n 

s A {t) = A A ( t ) + XX(') + E*f (*)> 

i=l j=l 
m n 

S B {t) = A B (t) + Xi B (t) + E Y *®> 

i=l 3=1 

where A (t) and A B (t) are seasonal floors, X A and Xf are common OU processes, i.e. they are 
driven by the same jump process Lj. A different approach based on copula is proposed in [5] where 
the joint evolution of electricity and gas prices is modeled by a bivariate non-Gaussian OU pure 
jump process with a non-symmetric copula. 

In this paper, we propose an alternative class (arithmetic and geometric) of models to reproduce 
adequately the statistical features of gas and electricity spot prices based on parameterized local 
volatility processes. The spiky behaviour of both spot prices is captured without introducing jump 
diffusion processes. More precisely, the deseasonalized (log) prices processes are modelized by the 
sum of an Ornstein-Uhlenbeck process (which is common for both commodities) and an independent 
stationary diffusion process. The construction of this stationary diffusion is similar to [7] where 
diffusion models with linear drift and prespecified marginal distribution are investigated with an 
application in a different context. The selected parameterized diffusion coefficient allows to capture 
"cluster of volatility" (e.g. period of high volatility which implies prices spikes). 
Moreover, this approach provides a significant advantage over the class of jump diffusion models 
since the calibration process involves only classical statistical tools like least squares method and 
maximum likellihood estimations so that it is robust and fast. It allows to reproduce (for the first 
time to our knowledge) both the auto-correlation and the cross-correlation strutures between two 
energies. The model was successfully tested on several markets and seems to fit well the statistical 
features and the marginal distributions of gas and electricity spot prices. 

Our results are presented as follows. Section 2 is devoted to the description of the stylised 
features of gas and electricity spot prices. Then, in Section 3, we briefly recall some important 
theoretical results on which are based our model. To be more precise, we recall how to construct 
a mean reverting diffusion process X solution of a stochastic differential equation (SDE) with a 
prespecified continuous invariant density /. Such diffusions involves parameterized local volatility 
processes. In Section 4, we present the model of our choice and focus on the calibration procedure. 
In the last section, we perform the calibration on the data sets coming from the NBP for the 
gas spot price and the Powernext market for the electricity spot price. Then, we proceed to the 
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simulation and, finally, analyze the impact of the modelization by measuring the risk of an energy 
related portfolio using several models. We show that introducing the cross-commodity correlation 
structure can greatly modify the risk of a portfolio. 



2 Stylised features of gas and electricity spot prices 
2.1 Seasonality 

A first characteristic of gas and electricity (and many commodities) prices is the presence of annual 
(and possibly multi-time scales) seasonality and a trend (see e.g. [12], [20]). For each commodity, 
we model the seasonality and the trend component of the logarithmic spot prices with the mean 
level functions around which spot prices fluctuate 

log g(t) = a 9 + b 9 t + ^ c l cos ( T~ ) + d k sin 

k=l 



Ik J V h 



log e(t) =a e + b e t + J2 4 cos (j^J + d% sin 



where Ik = |_252//cJ, k = 1, ...,m, [^J denotes the integer part of x. For instance, if we choose 
m = 2, we only consider a seasonal function over the year and the semester. We assume 252 trading 
days in a year except for electricity spot price on Powernext which has 365 trading days in a year 
so that, we have to take into account this particularity in the seasonality function. The coefficients 
above are estimated using ordinary least squares. The log-seasonality functions are represented with 
the estimated values for m = 2 using gas spot price at the NBP and electricity spot price from the 
Powernext market in Figure [TJ All parameters are not significant at the 5% level. We only report 
and take into account the significant values 1 . We checked the seasonality over week, month and 
quater, but the coefficients were not significant. 




Figure 1: The fitted log-seasonality functions log(g(i)) and log(e(t)) 



Now we focus our attention on the deseasonalized data Y 9 (t) := log S 9 (t) — log g(t) and Y e (t) := 
log S e (t) — loge(f) for the specification of the model. A geometric model consists in modelling 
the stochastic processes Y 9 (t) and Y e {t) whereas an arithmetic model consists in modelling the 
stochastic processes e Y3 ^ and e Ye ^\ 

V = 1.53,6 s = 0.000688, c? = 0.121, d\ = 0.0287, cf = 0.00533 et a e = 3.02, b e = 0.000405, cf = 0.138, dl = 
0.0368. 
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2.2 Spikes and heavy tails 



Electricity has very limited storage possibilities. It induces the possibility of spikes in spot prices. 
Natural gas can be stored but it is often costly, so that it shares the spiky behaviour of spot 
electricity prices. Gas and electricity markets share this similarity as it can be seen in Figure 
presenting the electricity spot prices coming from the Powernext market on the left and gas spot 
prices at the National Balancing Point (NBP) on the right. From a stochastic modelling point of 
view, spikes are commonly represented by jump diffusions with mean reversion. However (to the 
best of our knowledge) there is no evidence that it is rather jumps than spikes caused by clusters 
of volatility for instance. 




Figure 2: Electricity spot prices on the Powernext market (on the left) and gas spot prices at the 
NBP (on the right) for the period 14 January 2003 till 20 August 2008. 



The histograms of Y g and Y e with the fitted normal density curve is presented in Figure 
We observe that the two residuals time series Y 9 and Y e are far from being normally distributed. 
The excess of kurtosis of Y g and Y e are respectively equal to 4.5 and 2.3 meaning that the two 
distributions are peaked and have heavy tails. The skewness of Y 9 and Y e are respectively equal 
to 0.77 and 0.57 meaning that the two distributions are not symmetric. 



Histogram of Yg Histogram of Ye 




Figure 3: Histograms of Y 9 and Y e with normal density curves. 



2.3 Mean reversion and long term dependency 

Gas and Electricity spot prices are known to be stationary. This can be tested using an augmented 
Dickey-Fuller test (ADF) or the Phillips-Perron test. For the UKPX, Powernext electricity spot 
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prices and gas spot prices at the NBP the unit root hypothesis was rejected using both tests. 
Figure [5] shows that gas and electricity deseasonalized prices are strongly linked by a long term 
dependency, i.e. it seems that there is a stochastic equilibrium between Y 9 (t) and Y e (t) from 
which they cannot deviate for a long time. This long term dependency can be observed on the 
cross-correlation function. 




Figure 4: The log-deseasonalized gas (normal line) and electricity spot (dashed line) prices 



2.4 Auto-correlation and cross-correlation 

In energy spot price modelling, the auto-correlation functions (ACFs) are often analyzed. The ACFs 
of both Y 9 (t) (respectively e Y3 ^), p 9 , on one hand Y e (t) (respectively e Ye ^), p e , on the other hand 
present both a two-scale (or three-scale at most) decreasing behaviour with one quickly decreasing 
component and one or two slow decreasing components. The same behaviour is observed on the 
cross-correlation function (CCF) p 9 ' e . This kind of decreasing ACFs and CCF are well explained 
by sum of decreasing exponentials components, namely for r > 0: 

p 9 ( T ) = Corr (Y 9 (t + r), Y 9 (t)) = <^e~ A ? T + (1 
p e ( r ) = Corr (Y e (t + r),Y e (t)) = <f x e~^ T + (1 
^{t) = Corr {Y 9 {t + r), Y e {t)) = ^e^'^ . 

For the sake of simplicity in our stochastic modelization, we focused on one type of cross- 
correlation Corr (Y 9 (t + r), Y e (t)) and we assumed that the cross-correlation is symetric that is 
Corr (Y 9 (t + r), Y e (t)) = Corr (Y e (t + r), Y 9 (t)) which is a rather natural approximation. We 
observed that the slower rates of mean reversion for each commodities are quite similar A| = A| 
and that a rather good approximation is obtained by setting A 9 ' 6 = Af, = A^. Using a least squares 
approach, we fitted simulteanously p 9 {r), p e (r), p 9 ' e (r) (r = 1, 150) to the empirical ACFs and 
CCF. We assumed that the observed spot prices have reached the stationarity. Both empirical and 
fitted ACFs 2 and CCF 3 are plotted in Figure El 

We can see the separation into a fast speed of mean reversion for gas and electricity spot prices Af 
and Af which corresponds to a correlation dependence of approximately 2 and 30 days probably due 
to the spikes components whereas the slower speed of mean reversion corresponds to a correlation 

2 4>\ = 0.43, Xf = 7.2, and 4>t = 0.49, Af = 69.4 
3 4> 9 ' e = 0.53, Af = AS = A 9 ' e = 2.6 



-ft)*-*, 
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(a) ACF of Y 9 (b) ACF of Y e (c) CCF of (V 9 , Y e ) 

Figure 5: Empirical ACF and CCF of deseasonalized gas spot price and electricity spot price 



dependence of 64 days and corresponds to the stochastic equilibrium or the normal variation of gas 
and electricity spot prices. 



3 Theoretical background 

In order to modelize heavy tails (and spikes) of stationary spot prices distribution, a natural idea 
is to consider an ergodic diffusion process like representation of deseasonalized spot prices. 

In this section, we briefly recall how to construct a one dimensional process X solution of a 
stochastic differential equation with a prespecified continuous invariant density /. Throughout the 
sequel we assume that / is a strictly positive bounded continuous probability density on (l,r) (and 
zero outside (I, r)). 



3.1 The general case 

Let (Xt) t>0 the diffusion solution of the following stochastic differential equation (SDE) 

dX t = b(X t )dt + a(X t )dB t , X E (I, r), (E b>a ) 

where b : (l,r) — > R and a : (I, r) — > R are locally Lipschitz functions and that a is not degenerate 
on (l,r) i.e. \/x E (l,r),a 2 (x) > 0. We introduce for the diffusion (Xj) i>0 , the scale function 
p : (I, r) — >■ R defined for any c S (/, r) by 

f x ( f y 2b(z) \ 
Vx G (/, r), p(x) = J exp \ - J ^2^y dz J d V> 

and the speed measure density m : (l,r) — > R!j_ defined by 

Vx e (l,r), m (x) = p/(x) l 2{x) = ^yexp |||d.) . (1) 



We recall (see e.g. [HI [T7]) that the process \^p(X^)j with ( = inf {t ^ 0, X t = I or X t = r] is 

a local martingale if and only if p is the scale function (unique up to an affine transformation). 
Moreover, if the diffusion (Xt) t>0 is positive recurrent, the stationary probability distribution v 
defined on (/, r) satisfies 

v{dx) = Cm(x)dx with C = (^j m{x)dx 

This classical result is the key to construct a one-dimensional ergodic process that fits prescribed 
stationary probability distribution. For a more general result to construct an inhomogeneous Markov 
martingale process that has prespecified marginal density we refer to [19|. 
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Proposition 3.1. Let b : (I, r) — > R be a continuous drift function. Suppose that b and f satisfy 
the following conditions 



Vx€(l,r), / b(y)f(y)dy>0, and / b(y)f(y)dy = 0, (H B ) 
Ji Ji 

Then there exists a unique continuous diffusion function defined by 



Vx G (I, r), a(x) 



J l x b(y)f(y)dy 



/(*) 



such that {Eb i(T ) has a unique solution {Xt) t>0 , which is an ergodic diffusion process with stationary 
distribution v satisfying v{dx) = f(x)dx. 

Further details of the proof outlined below can be found in [7J. 

Proof. Let B be the function defined by B(x) = b{y)f{y)dy. One checks easily that the scale 
function of [Xt) t>(j satisfies 



f x 1 

VxG(/,r), p(x)=B{c)J ^y d 2/- 



One then obtains that linXj.^ p[x) = — oo and lim x _). r p[x) = +oo. 

On the other hand, the speed measure of (X t ) t>Q has density m that satisfies 

Vx G (t, r), m[x) - 



B(x)p'(x) B(c)' 

The normalized speed measure density is then equal to the probability density /. 

To prove existence and uniqueness of the solution (-Xt) t>0 > one proves existence and uniqueness 
of the process (p(Xt)) t>0 satisfying a SDE without drift (see [E]). □ 

Corollary 3.2. Let b : x G (I, r) \— > —X(x — fi) and assume that probability density f has ex- 
pectation fj, and finite variance. Then there exists a unique continuous diffusion function defined 
by 



Vx G (I, r), a(x) 



t f l x 2\(y-y)f(y)dy 



such that [E^ a ) has a unique solution {Xt) t>Q , which is an ergodic diffusion process with stationary 
distribution v satisfying v(dx) = f(x)dx, and ACF given by 

Vt,r^0, coic(X t+T ,X t ) = e~ Xr . 

The squared diffusion coefficients are explicitly known for a large number of commonly used 
probability diffusions. However, for some specific distributions, it is not possible to obtain a closed 
form of the diffusion coefficient. An approximation based on saddlepoint technique and the moment 
generating function (which is generaly known explicitly) is developed in [7]. 



3.2 Quasi-Saddlepoint approximation 

We first recall that saddlepoint approximations are constructed by performing various operations on 
the moment generating function (MGF) of a random variable (see e.g. [8]). Let X be an absolutely 
continuous random variable with density / (with respect to the Lebesgue measure on (l,r)), moment 
generating function M(t) and cumulant-generating function n(t) = logM(i). Then the first-order 
saddlepoint density approximation to / is given by 

Vx G (l,r), f{x) = {27TK"(i x )y 1/2 e-( { ^- K ^), 
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where t = t x is the (unique) solution to the saddlepoint equation n'(t) = x, and primes denote 
derivatives. We assume that the probability density / has expectation fj,, i.e. jjL= k'(0). 

Considering the continuous differentiable function t : x i— > t x , an integration by parts gives 



t(y)dy = i(x)x - / i'(y)ydy, 
Jo 

rx 

= i{x)x - I dn(i(y)), 
Jo 

since y = K r (t(y)). The saddlepoint density / writes then 

Vx G (l,r), /(x) = (2W(£(x)))- 1/2 exp (- £ t(y)dy) . (2) 

To construct an ergodic process (X t ) t>(j solution of (E^ a ) with prespicified stationary density /, 

the exponential terms that appear in ([2]) and ([1]) suggest the relation = t. This construction is 
not exact but in [7J is proved that the speed density m of X is approximately propotional to the 

saddlepoint density /. To be precise both \/ n"(i(x)) and a 2 (x) are approximately proportional to 



k"(0) + (O)t(x) near the mean of the distribution. From now this normalized speed density m 
will be called the quasi- saddlepoint density approximation to /. 

To summarize, if the saddlepoint function t is explicity known and efficiently computed, then we 
consider the diffusion with drift b, such that b > on and b < on (/x,r), and with diffusion 
coefficient 



Vx G (I, r), a(x) 



-26(g) 



which is ergodic with stationary distribution fix) = e K (*( x ))) ( w here c is a normalizing 

factor), the quasi-saddlepoint density approximation to / (see [7] Theorem 3.1 for more details). 

The following example will become useful later when we are going to modelize deseasonalized 
gas and electricity spot prices. 

Example 3.1. The NIG- distribution The normal-inverse Gaussian (NIG) distribution is a member 
of the class of generalized hyperbolic distributions (see e.g. [3]). The NIG density is given by 



aSKi aJ5 2 + (x - I) 2 . 

f{x) = V / I x rSy/^+Kx-l) x£R 

TT^d 2 + (X - iy 

where /3 G R, a > |/3|, 5 > 0, / £ R and K\ is the the modified Bessel function of third order and 
index 1. Note that if X ~ NIG (a, (5,8,1) then its two first moments are 

E[X] = / + ^=^= et var(X) = — - — - T . 



V 7 " 2 - /3 2 (a 2 -/? 2 ) 
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The two parameters 5 and Z determine respectivelly the scale and the location of the law, and the 
two parameters a and f3 determine the shape: a being responsible for the tail heavyness and /3 for 
the skewness (asymmetry). 

The cumulant-generating function k of the NIG distribution is defined for all t such that \/3+t\ < 
a by 

K (t) = lt + 5 [^a 2 - /3 2 - y/a 2 - (/? + t) 2 ) , 
and the saddlepoint function is defined by 

Vx G R, t(x) = a ( x ~ 1 ^ _ f3. 

s/p + (x - iy 



8 



In order to have an Ornstein process solution of (E^ a ) with stationary density the quasi-saddlepoint 



density approximation / to /, we consider the following drift and diffusion functions 



VxGR, b(x) = -X(x-fJL) and a\x) = + {x , Z {x =L, (3) 

4 Cross-commodity multi-factor model 

In this section, we present two class of cross-commodity multi factor models: the geometric and 
the arithmetic class. Those two class are commonly used in stochastic modelling of commodity 
prices. The first one ensures the positivity of simulated spot prices. However, when dealing with 
forward contracts which have a delivery period or options pricing, the second one is analytically 
more tractable. Both class of models are based on stationary diffusion-type models analyzed in [7] . 

4.1 Proposed modelization 

In this section, we propose an alternative model which captures the stylized features described in 
Section 2 without introducing jump diffusions. Sums of this kind of diffusions can fit the multi-scale 
ACFs and the CCF obtained for the deseasonalized gas and electricity spot prices. 

In order to represent the ACFs and CCF of gas and electricity deseasonalized spot prices, we 
are led to introduce stochastic processes that are sums of diffusions defined by (E^). To be more 
precise, we focus on the following two factor modelization for the deseasonalized log spot prices Y 9 
and Y e 

Yf = Xf + Z t , and Y t e = Xf + Z t , (4) 
where (Zf) f>Q , (Xf) t>Q and (Xf) i>Q are mutually independant processes defined as following: 

• the process (Zt) t>0 accounts for the stochastic equilibrium between both commodities with a 
slow rate of mean reversion X z = A2 = A|- Thus, it represents the normal variation and will 
be defined by an Ornstein-Uhlenbeck process 

dZ t = -X z Z t dt + a z dW t z , (5) 

with X z > and a z € R. Note that Z is ergodic with the Gaussian invariant probability 
Af(0,cr 2 j2X z ). 

• the processes (Xf) 4>Q and (Xf) i>Q represent the spikes component for each commodity. We 
modelize them by general Ornstein processes with high rate of mean reversion X g = Af > 
and A e = A| > 0, namely 

dX> = -Xj (Xi - H ) dt + aj(xi; 0,-)dW/ 5 j = g, e, (6) 



where Oj is a parametric diffusion function such that {Xf) t>Q is an ergodic diffusion with 
invariant probability f J (.,0j). 

Remark 4.1. The following contraction can be extended to a more general multi- factor model. 
We can consider m general Ornstein processes and p Ornstein-Uhlenbeck processes so that 



Y9{t)=Y J m)+i2 z j^ 

i=i 3=1 

rn p 

Y*( t ) = j2m)+T, z i( t )> 

i=i j=i 

where all processes are assumed to be mutually independent, i.e. driven by independent Wiener 
processes. We already observed that a two-factor model (m = 1 and p = 1) fits the ACFs and CCF 
well. 
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Proposition 4.1 (The correlation structures). Let Y 9 , Y e be the processes defined in Then, 
the A CFs of Y 9 and Y e with lag r > are given by 



p 9 (r) = cor (Yf +T ,Yf) = 4> g e- x ^ + (i _ ^)e~^\ 

p e (r) = cor {Yf +T ,Yf) = <p e e~ x ^ + (1 - e )e~ A ^, 

where 

_Xax{X 9 (t)) _Var(X e (i)) 

09 " Var(F^))' an(i ^" Vax(y«(t))' 

77ie CCF i/rai/i Zag t > is given by 



p* e (r) :=cor (F/ +T ,y t e )- 



A /Var(y9(t))Var(y 



(*)) 



From the definition of d> ae . we find that o\ = 2A 2 (/> 9 , eX /Var (Y 9 (t)) Var (T e (f)), where the last 
term is the product of the two stationary variance of the two processes. Consequently, one can 
easily derive o z from the ACFs and CCF calibration. 



4.2 Calibration 

We propose a three-step calibration procedure for the model described above. 



Step 1: Deseasonalizing spot prices 

We fit the seasonality functions g(t) and e(t) defined in section 2.1 to the logarithmic spot prices. 
The parameters of the functions are estimated using the least squares approach. Now, we focus on 
the deseasonalized spot prices Y 9 and Y e defined by 

Y 9 (t) = \og (S 9 (t))- log (g(t)) and Y e (t) = log (S e (t)) - log ( e (t)) . 

One can consider the deseasonalized spot prices e Y9 ^ and e Ye ^ instead of this geometric approach. 



Step 2: ACFs and CCF 

The least squares method consists in fitting the empirical ACFs p 9 (r), p e {r) and CCF p 9,e (r) defined 
in section 2.4 to the empirical ones (p 9 (r)) r=1 ; , (/5 e (r)) r=1 t , (p 9 ' e (T)) T=1 l in order to derive 
the three speeds of mean reversion Af, Af, A^ with the diffusion coefficient o z of the stochastic 
equilibrium process Z . This can be done by minimizing the sum of squared differences, namely 

i 

argmin £ ((p 9 (r) - ^(r)) 2 + (p e (r) - p e (r)) 2 + (^ e (r) - p^{r)f) . 

X g ,\e,X z ,a z j 

Stability tests showed that the estimates are robust with respect to small changes in the initial 
values of the parameters. 



Step 3: Estimating the parameters of the spikes component 

The final step consists in statiscally estimating the parameters 6 g of the invariant density f 9 (.,9 g ) 
of the process X 9 and the parameters 9 e of the invariant density f e (.,8 e ) of the process X e . For 
instance, if one decide to choose the quasi-saddlepoint approximation to the NIG density for f 9 and 
/ e , there will be four parameters to fit for each density. The model proposed is a sum of diffusion 
processes and hence is not Markovian. Thus, the likelihood cannot be written down explicitly. To 
overcome this problem, we use the maximum likelihood estimation of order m (m = or m = 1 in 
our case) method for stationary processes introduced in [1] . Strong consistency and a Central Limit 
Theorem are proved for such estimates. It consists in approximating the log-likelihood of the serie 
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(yi)i<fc<n (j = 5) e )> where n is the number of sample points, by a sum whose generic term is the 

4 



density function of Yl conditional on the m most recent observations, for some m > 0, namely, 



= £ log (^Vjj/fH) 



(7) 



fc=i 



where y£ m := (y k _ m ,--- illl-i) an d | Vk™',®]) is the conditional probability density function 
of Yl given Y£' = y k m for the parameters 6j. Note that if m = 0, there is no conditioning 
and ^ is simply the marginal density of Y£, k = 1, • • • , n, which is the convolution of and 
X^. We suppose that (^)i<fc<n (resp. {Zk)\<k<n) is ergodic with stationary distribution f J (.;9j) 
(resp. with Gaussian invariant probability M (0,a~z 2 ■= o"^/2A z )) so that the conditional probability 
density function is given by 



W fc ;*i) 



7T 



d«, j = g, e. 



(8) 



Note that it corresponds to the case of (^?)i<fc< n is independent and identically distributed random 



variables having the distribution of the stationary distribution of Y 3 . 

Numerically, the above integral can be approximated using a Gauss-Hermite quadrature method, 
namely 

n / 



1 n 



k=l 



Uk',Vj }w k , j = g,e, 



where (uk)\^k^n are the roots of the Hermite polynomial P n and {wk)i <k<n 
weights given by 



are the associated 



k = 1, n. 



If m = 1, we need to compute the transition probability density p y j ^ Y j_ y) _(.; 0j) of (Y/) t>0 for 

j = g,e. The two series (^)i<fc<n an d {Zk)\<k<n ar e discrete observations of ([6]) and ©. Let g 
be a Borel bounded test function and k £ 1, ■ ■ ■ ,n — 1, by conditioning we have 



E 



9(Y£ +1 ) I Yl = y k 



/ 



9( Y Li) I Y k = Vk, Zk = z P [Z k = z | = y) dz 



g{x{ +1 + Zfc+i) | = y k , z k 



z k = z | y, J = y ) d^ 



Note that the two processes X J and Z are independent so that if we denote by p x 
Px(tk, tk+i, x k , .) and pz k (z k , .) := pz(t k ,t k+ i, z k , .), the conditional probability density functions of 
X J k+1 and given Z k = z, X{ = x J k , the expectation E g(X J k+l + Z k+l ) \ Y k = y J k , Z k = z 

given by 



is 



/ g( u ) / p x j(yl~ z,v)p Zk (z,u- v)dvdu. 



Moreover, we have 



Z k = z I Y£ 



k ~ Uk 



(Z k = z, Yl = y{) _F(Z k = z)F(x{ = y{ - z) 



j 

Vk 



{yl = A) 



where, 



Yl 



k — y°k 



+oo 



2'ndz 



-e 25 z du, 



11 



and P ( X 3 k =yi> k -z\ = ft (y{ - z; 6j J , P (Z k = z) = ^=^e ^ . Finally, one easily identifies 
the transition probability density p y j ,y,_ j(y; Oj) which is given by 

1 ft(y{-z;0 3 )e~^ z2 

Px j (yi ~ z > v ^ z k. ( z > u - v ) dv -t^z r~, — ?\ du - 



1 »2 



Note that we have pz k (zk, z) = e 2rT z , with a z = a z\J 1 e 2 \ z Z using an exact scheme of 

the Ornstein-Uhlenbeck process (Zfc)i<fc<n of step A > 0, namely 



1 _ e -2A z A 

Z k+1 = e- x * A Z k + a z J 2 ^ G* k+1 , (9) 



where {G k ) k>1 is a sequence of i.i.d. standard normal random variables. However, in most cases, 



there is no closed expression for p x j(x k ,.). To overcome this problem one solution is to consider 



the transition probability density function p x j(x k , .) of the Euler scheme ( X k 



k>0 



X{ +1 = e~^X{ + H (l - e^ A ) + (Tj (Xl; 6 3 ) ^ 1 ^ G{ +v k > (10) 
where (Gk*) ^ s a sec iuence of i.i.d. standard normal random variables independent of {G k ) k>v 

1 „2 



SO 



Consequently, Pj^K,.) = ^.^^.f "V ^ , with aj(x 3 k ;9j) = aj (x> k ,0j 
that we have 

p x Ay{ - z,v)p Zk (z,u- v)dv = —— e^vl^ {U - mk} , 

k V2Tra(y J k ,z;03) 

where for k £ {1, • • • ,n}, m k = e~ XzA z + e~ x ^ A (y 3 k — z) + fij(l - e~ x ^ A ) and a 2 (y{, z;6 j ) = 

Remark 4.2. In [10], a transition probability density function based on Milstein scheme is used. 
In [18], a gaussian transition probability density function with Taylor expansions is used to propose 
an efficient estimator for Oj. 

The method of maximum likelihood of order m estimates Qj, m by finding the value of 6j that 
maximizes ([7]) using standard numerical optimization procedure. 

5 Simulation and application 

5.1 Empirical results on Powernext and NBP spot prices 

In this section, we perform the calibration procedure on electricity spot prices coming from the 
Powernext market and on gas spot prices at the NBP. Then, we perform a simulation with estimated 
parameters over the same period. To avoid negative prices, we choose to represent spot prices by 
an arithmetic model, namely 

S°(t)=g(t) xe X9 W +z W, (11) 
S e (t) = e(t)xe xe ^ +Z ^\ (12) 

where g(t), e(t) are the trend and seasonality functions defined in Section \2. 11 X 9 , X e are solutions 



of (Ef, >CT \ with b and a defined in ([3]) and Z is a Gaussian Ornstein-Uhlenbeck process solution of 
©■ 
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We choose the NIG distribution for those two processes in order to capture the heavy tails 
behavior observed on data, i.e. large values with low probability that cannot be obtained by a 
Gaussian process. We observed that the quasi-saddlepoint approximation of the NIG-distribution 
is well suited to represent the two spike components. One can choose another distribution and devise 
the same calibration process as in the previous section. The results of steps 1 and 2 of the calibration 
procedure are reported in Figure [1] and the quality of the ACFs and CCF fits is represented in Figure 
[SJ Now, we proceed to the estimation of the four parameters 6 g = (a g , /3 g ,5 g ,l g ) of the process 
X 9 and the four parameters 9 e = (a e , f3 e ,5 e ,l e ) of the process X e using the maximum likelihood 
estimation method described in the previous section on the deseasonalized spot prices. We observed 
that the maximum likelihood estimation method of order 4 is more robust and gives better results 
than the one of order l 5 . The initial parameters are set to (1, 0, 1, 0) for both components. 

The algorithms converged quickly. The diffusion coefficient functions (fj(., 9j), j = g, e, with the 
fitted parameters, are documented in Figure [6j We see that the shape of the diffuion coefficients 
are quite similar for the gas and electricity spot deseasonalized spot prices. Spikes are obtained 
when the processes Y 9 and Y e are far from their mean by clusters of volatility, i.e. periods of 
high volatility. As we see, large values are more likely and the asymmetry is more pronounced for 
electricity spot prices than for gas spot prices. We clearly see spikes as cluster of volatility are more 
probable and more intense for electricity deseasonalized spot prices than for gas deseasonalized spot 
prices. 



Diffusion coefficient for gas Diffusion coefficient for electricity 




Figure 6: Squared diffusion coefficients using fitted parameters with maximum likelihood estimation 
of order (normal lines) and of order 1 (dashed lines). 

In order to simulate price trajectories, we consider the Euler-Maruyama schemes defined by 
(|10p and ([9]). If one is concerned by estimating some quantities (for instance quantiles) on only 
one trajectory then one should replace the above Euler schemes of X 9 and X e with their respective 
Milstein schemes X 9 and X e in order to achieve a smaller strong error rate. It consists in devising 
the following schemes for j = g,e, 




where a/ is the first derivative of Uj. 



4 Fitted parameters of order are: a g = 1.93, f} g = 0.90, S g = 2.25e - 3, l g = -8.8e - 3 and a e = 3.49, f} e = 1.24, 

S e = 0.08, l e = 0.11. 

5 Fitted parameters of order 1 are: a g = 0.76, /3 g = 7.8e — 2,S g = 7.8e - 4, l g = -0.11 and a e = 1.56, p e = 0.34, 
S e = Lie -2, l e = 0.16 
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In the following simulations, we consider Milstein schemes of step tk = kA, with A = 
Next, we add to the simulated processes the two seasonality functions. In Figure [71 the simulated 
deseasonalized spot prices are represented. We see that both commodities are strongly linked and 
that the model mimics the statistical behaviour of the deseasonalized spot prices. In Figure El the 
simulated spot prices are represented. In Figure El both simulated and historical ACFs and CCF 
are plotted. We clearly see that the model reproduces the correlation structures. 




Figure 7: A simulation of gas (normal line) and electricity (dotted line) deseasonalized spot prices. 




Jan 03 Jan 04 Jan 05 Jan 06 Jan 07 Jan 08 Jan 03 Jan 04 Jan 05 Jan 06 Jan 07 Jan 08 



Figure 8: Simulated Electricity spot prices on the Powernext market on the left and Gas spot prices 
at the NBP on the right for the period 14 January 2003 till 20 August 2008. 



5.2 Application: measuring risk of a cross-commodity portfolio 

In this section, we aim at measuring the risk of a portfolio composed of a short position in a power 
plant that produces electricity from gas day by day t± < ... < t n for several maturities T = tjy = 6 
months, 1 year and 3 years. The loss at time of the portfolio with a time horizon T can be written 

N 

L T = Y j e~^ (St k -h R S? k -C) + -Pt, 

k=l 

where r = 5% is the annual interest rate, = 3 denotes the Heat Rate, C = 5 €/MWh denotes 
the generation costs and where P£ is an estimation of the price of the option on the power plant 
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(a) ACF of one simulation of Y 9 



(b) ACF of one simulation of Y e 



(c) CCF of the simulations 

Figure 9: ACFs and CCF of simulated gas and electricity spot prices (normal lines) with the 
historical ACFs and CCF (dotted lines). 



obtained by a crude Monte Carlo simulation, namely 



pc 



N 

E' 

k=i 



-it 



fc E 



(Sf k -h R S 9 tk -C)_ 



Since gas and electricity markets are incomplete, we price and estimate risk measures under the 
historical probability. In order to measure the risk, we consider the Value-at-Risk (VaR), which is 
certainly the most commonly used risk measures in the context of risk management. By definition, 
the Value-at-Risk at level a 6 (0, 1) (VaR a ) of a given portfolio is the lowest amount not exceeded 
by its loss with probability a. In this example, we set a = 95%. Actually, for the considered 
portfolio, the VaR a is the unique solution £ of the equation 



P [Lt < £] 



O'. 



The portfolio's VaR a is just a quantile of its loss and is interpreted as a reasonable worst case level. 

Now, we are interested in measuring the impact of the proposed model for gas and electricity 
spot prices on the portfolio's VaR. In order to do that, we consider three different models: 

• Case 1: the mean-reverting cross-commodity model (in its geometric form) proposed in this 
paper and defined by ([TT]) and (fT2j) . It modelizes typical features of gas and electricity spot 
prices likes spikes and the long term dependency. 
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Maturity 


pc 


(±Error) 


VaR a 


Case 1 

iriOpOSCQ IIKJUei 1 


6 months 


83.3 


(±3.3) 


262.4 


1 year 
3 years 


220.1 
745.0 


(±5.5) 
(±11.2) 


495.4 
1081.0 


Case 2 

(No cross-correlation) 


6 months 


51.2 


(±2.9) 


250.1 


1 year 
3 years 


222.6 
850.6 


(±8.4) 
(±21.3) 


880.2 
2213.1 


Case 3 

(Gaussian model) 


6 months 


32.9 


(±1.1) 


107.7 


1 year 


129.8 


(±2.7) 


275.9 


3 years 


437.1 


(±5.8) 


565.5 



Table 1: Estimation of the price of the Power plant and the VaR a of the portfolio. 



Case 2: a slight modification of the previous model in which we do not take into account 
the dependence of the two energy spot prices. To be more precise, we consider the following 
model specification 

S 9 (t)=g(t) xe^W+^W, 
S e (t) = e(t)xe xe V +ze V, 



where X 9 and X e are solutions of (Eb^) with b and a defined in ([3|), and where Z 9 , Z e are 
two independant Gaussian OU processes solution of ([5]). By this model, we want to measure 
the impact on the VaR a of the long term dependency modeling. The calibration process is 
slightly modified since S 9 and S e are now independent. The step 2 is replaced by two different 
minimizations corresponding to each ACF. Steps 1 and 3 remain unchanged. 

• Case 3: a slight modification of the case 1 in which we do not modelize the spikes feature. To 
be more precise, we replace the NIG-distributed processes by Gaussian Ornstein-Uhlenbeck 
processes, namely 

S 9 (t)=g(t)xe Z3 ^ +z( - t \ 
S e (t) = e(t)xe z ^ +Z ^, 

where Z 9 , Z e , Z are three different Gaussian OU processes solution of ([5]). By this model, 
we want to quantify the impact on the VaR a of the spike feature of gas and electricity spot 
prices. 

In each case, we estimate Pj, and the VaR a using 10 000 Monte Carlo simulations. We devise 
Euler schemes of step t k = kA with A = ^ • I n order to estimate the VaR Q , we use the inversion 
of the simulated empirical distribution function. 

Remark 5.1. Since gas and electricity spot prices are sums of diffusion processes solution of (E^,^), 
one can easily use the method investigated in [11] to estimate the VaR a and other risk measures. It 
is based on stochastic approximation algorithms with an adaptive variance reduction tool (uncon- 
strained importance sampling algorithm). The method is known to achieve good variance reduction 
when a ~ 1 as it is often the case. For the sake of simplicity, we only considered the classical 
method based on the inversion of the empirical distribution function. 

The results are summarized in Tables [TJ Note that for each case, the estimations are computed 
using the same pseudo-random number generator initialized with the same seed. The number in 
parentheses refers to the 95% confidence level. 

We observe that there are slight differences in terms of the price P£ between the case 1 and 2 
but huge differences in terms of risk. Taking into account the long term correlation between gas and 
electricity spot prices can reduce substantially the risk of this portfolio. Modeling independently 
each energy spot prices leads to an overestimation of the VaR a of the portfolio's loss. The results 
obtained by using the model investigated in case 3 shows that introducing the spikes behavior into 
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the model can increase greatly both Pj, and the risk of the portfolio. We also estimated the same 
quantities using the arithmetic version of the three models presented above. We obviously obtained 
different values from the ones presented but the same conclusions hold: modeling adequatly the cross 
correlation between gas and electricity spot prices reduces the risk of portfolio whereas modeling 
adequatly the spiky behavior of both commodities increases greatly the price of the option and the 
risk associated to the portfolio. 
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